#06/10/2012 Tai-Hsien Ou Yang


args=(commandArgs(TRUE))
cat(args,"\n")

if(length(args)==0){
    print("No arguments supplied.")
    ##supply default values
    k = 1
}else{
    for(i in 1:length(args)){
         eval(parse(text=args[[i]]))
    }
}

cat(k,"\n")

cat("=========================================\n")
cat("=====Chromosome-wise IBD Consensus=====\n")
cat("=========================================\n")
cat(k,"= Chromosome\n")
cat(getwd(),"\n")

#dataset<-scan(FILE_INPUT , list(0,"",0,"",0,0,0,"","",0,0,"",0,0,0))
#    dataset[["fid.1"]]=dataset[[1]]
#    dataset[["id.1"]]=dataset[[2]]
#    dataset[["fid.2"]]=dataset[[3]]
#    dataset[["id.2"]]=dataset[[4]]
#    dataset[["chr"]]=dataset[[5]]
#    dataset[["seg.s.bp"]]=dataset[[6]]
#    dataset[["seg.e.bp"]]=dataset[[7]]
#    dataset[["seg.s.snp"]]=dataset[[8]]
#    dataset[["seg.e.snp"]]=dataset[[9]]
#    dataset[["total.snp"]]=dataset[[10]]
#    dataset[["total.gene"]]=dataset[[11]]
#    dataset[["unit"]]=dataset[[12]]
#    dataset[["mismatch"]]=dataset[[13]]
#    dataset[["hom.1"]]=dataset[[14]]
#    dataset[["hom.2"]]=dataset[[15]]




FILE_LOCATION   ="/ifs/scratch/c2b2/ip_lab/to2232/"
FILE_INPUT_PED = "parchild.lst"

FILE_INPUT_IBD        = paste("Soapv2.",k,".match", sep="")
FILE_OUT_P            = paste("p.",k,".txt", sep="")
FILE_OUT_N            = paste("n.",k,".txt", sep="")
FILE_OUT_PR           = paste("p.",k,".Rdata", sep="")
FILE_OUT_NR           = paste("n.",k,".Rdata", sep="")



  setwd(FILE_LOCATION)

#Read Reference List
	 
 dataset.ped<-scan(FILE_INPUT_PED , list(ch="",pr=""))

#Read Raw Data  

  #dataset.pair<-scan(FILE_INPUT_IBD , list(id="",len.avg=0,len.var=0,len.total=0,ibd1=0,ibd2r=0))
  #dataset.pair<-scan(FILE_INPUT_IBD, list(fid=0,id1="",fid2=0,id2="",chr=0,segs.bp=0,sege.bp=0,segs.snp="",sege.snp="",total.snp=0,total.gene=0,unit="",mm=0,hom1=0,hom2=0))

  dataset.id1<-do.call("cbind", scan(FILE_INPUT_IBD, what = list(NULL,"",NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL), flush = TRUE))
  dataset.id2<-do.call("cbind", scan(FILE_INPUT_IBD, what = list(NULL,NULL,NULL,"",NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL), flush = TRUE))
  dataset.segs.bp<-do.call("cbind", scan(FILE_INPUT_IBD, what = list(NULL,NULL,NULL,NULL,NULL,double(0),NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL), flush = TRUE))
  dataset.sege.bp<-do.call("cbind", scan(FILE_INPUT_IBD, what = list(NULL,NULL,NULL,NULL,NULL,NULL,double(0),NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL), flush = TRUE))
  dataset.total.gene<-do.call("cbind", scan(FILE_INPUT_IBD, what = list(NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,double(0),NULL,NULL,NULL,NULL), flush = TRUE))
  dataset.hom1<-do.call("cbind", scan(FILE_INPUT_IBD, what = list(NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,double(0),NULL), flush = TRUE))
  dataset.hom2<-do.call("cbind", scan(FILE_INPUT_IBD, what = list(NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,double(0)), flush = TRUE))
  cat("Chr",k," loaded\n")

#Take out only what we want

parchild<-c()


    for(j in 1:(length(dataset.ped$ch)))
    {

	candidate.ch.f<-grep(dataset.ped$ch[j],dataset.id1)
	candidate.pr.f<-grep(dataset.ped$pr[j],dataset.id2[candidate.ch.f])
	if(length(candidate.pr.f)!=0)
	{	

		i=candidate.ch.f[candidate.pr.f[1]]
		pair<-cbind(dataset.ped$ch[j],dataset.ped$pr[j],dataset.segs.bp[i],dataset.sege.bp[i],dataset.total.gene[i],dataset.hom1[i],dataset.hom2[i],k,1)
	
	       parchild<-rbind(parchild,pair)
          	write.table(pair, FILE_OUT_P, row.names=F,col.names=F,append=T)

  		dataset.id1=dataset.id1[-i]
  		dataset.id2=dataset.id2[-i]
  		dataset.segs.bp=dataset.segs.bp[-i]
 		dataset.sege.bp=dataset.sege.bp[-i]
  		dataset.total.gene=dataset.total.gene[-i]
  		dataset.hom1=dataset.hom1[-i]
  		dataset.hom2=dataset.hom2[-i]


		cat("#",j,dataset.ped$ch[j]," ",dataset.ped$pr[j]," detected\n" )
	
	}

	#Also consider reversed condition
	candidate.ch.r<-grep(dataset.ped$ch[j],dataset.id2)
	candidate.pr.r<-grep(dataset.ped$pr[j],dataset.id1[candidate.ch.r])
	if(length(candidate.pr.r)!=0)
	{
	
		i=candidate.ch.r[candidate.pr.r[1]]
		pair<-cbind(dataset.ped$ch[j],dataset.ped$pr[j],dataset.segs.bp[i],dataset.sege.bp[i],dataset.total.gene[i],dataset.hom1[i],dataset.hom2[i],k,1)
	
	       parchild<-rbind(parchild,pair)
          	write.table(pair, FILE_OUT_P, row.names=F,col.names=F,append=T)

  		dataset.id1=dataset.id1[-i]
  		dataset.id2=dataset.id2[-i]
  		dataset.segs.bp=dataset.segs.bp[-i]
 		dataset.sege.bp=dataset.sege.bp[-i]
  		dataset.total.gene=dataset.total.gene[-i]
  		dataset.hom1=dataset.hom1[-i]
  		dataset.hom2=dataset.hom2[-i]

		cat("#",j,dataset.ped$ch[j]," ",dataset.ped$pr[j]," detected\n" )
		
	}

  
	}

	#Rfile
	save(parchild,file=FILE_OUT_PR)
	save(dataset.id1,dataset.id2,dataset.segs.bp,dataset.sege.bp,dataset.total.gene,file=FILE_OUT_NR)
